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We present an accurate ab-initio tight-binding hamiltonian for the transition-metal dichalco¬ 
genides, M 0 S 2 , MoSe 2 , WS 2 , WSe 2 , with a minimal basis (the d orbitals for the metal atoms and p 
orbitals for the chalcogen atoms) based on a transformation of the Kohn-Sham density function the¬ 
ory (DFT) hamiltonian to a basis of maximally localized Wannier functions (MLWF). The truncated 
tight-binding hamiltonian (TBH), with only on-site, first and partial second neighbor interactions, 
including spin-orbit coupling, provides a simple physical picture and the symmetry of the main 
band-structure features. Interlayer interactions between adjacent layers are modeled by transferable 
hopping terms between the chalcogen p orbitals. The full-range tight-binding hamiltonian (FTBH) 
can be reduced to hybrid-orbital k ■ p effective hamiltonians near the band extrema that captures 
important low-energy excitations. These ab-initio hamiltonians can serve as the starting point for 
applications to interacting many-body physics including optical transitions and Berry curvature of 
bands, of which we give some examples. 


I. INTRODUCTION 

The successful isolation of a single atomic layer of 
graphene 1 ushered in a new era in the study of two- 
dimensional materials, but its gapless band structure 
has hindered its applications in electronic devices that 
depend on the presence of a band gap. In contrast 
to graphene, layered transition metal dichalcogenides 
(TMDCs) 2-5 with chemical formula MX 2 (representative 
examples being M=Mo,W and X=S,Se), are semiconduc¬ 
tors with direct or indirect band-gaps that depend on the 
number of layers and are in the range of visible light (1- 
2 eV). These materials also display a whole spectrum 
of phenomena such as superconductivity 6 , magnetism', 
charge density waves (in TaS 2 ) observed by experiment 8,9 
and topological insulator phases predicted by theory 10,11 . 
For monolayer TMDCs, the broken inversion symmetry 
enables optical control of valley degrees of freedom 12 
and the strong spin-orbit coupling leads to valley-spin 
coupling 13 . The two valleys of the fc-points labeled K± 
can be manipulated by breaking the time-reversal sym¬ 
metry with the optical Stark effect 14 , by an external mag¬ 
netic field 15 18 or by a magnetic substrate 19 . Van der 
Waals heterostructures 2,20 with TMDCs or other two- 
dimensional layered materials such as hBN (an insulator) 
can also be fabricated. In heterostructures with incom¬ 
mensurate lattice periodicity, the induced interlayer po¬ 
tential from Moire patterns provides an additional way 
to control electronic properties 21,22 . With recent ad¬ 
vances in the experimental techniques for synthesizing 
these materials, novel optoelectronic, valleytronic and 
spintronic 13 applications have been proposed, that take 
advantage of the interplay between different features and 
the spin, valley and layer degrees of freedom. 

To address the whole range of interesting phenomena 
in the TMDC systems, it is crucial to have a simple and 
accurate model as a guide to the physics and the sym¬ 


metries involved. Thus far, only two types of approaches 
have been considered toward this goal: i) oversimplified 
models such as k • p theory and tight-binding hamiltoni¬ 
ans with few bands involving only the d orbitals of the 
metal atoms 13, 23-25 ; ii) fitted models that contain all the 
details of the ab-initio calculations, often parametrized 
as a tight-binding hamiltonian with a large number of 
parameters fitted to reproduce the ab-initio results 26,27 . 
While both types of approaches are valuable and interest¬ 
ing, they each have certain deficiencies: In the first type 
of approach, the physics is transparent due to the small 
number of basis functions and parameters involved, but 
not necessarily accurate or reliable. In the second type of 
approach, the fitted models have the required precision, 
but lack a simple basis and transparency in the physics; 
the large number of parameters involved in the fitting, 
which is not unique, may also lead to inconsistent results 
especially in relation to surface states as explained below. 

In this work, we present simplified models for TMDCs 
based on the full-range eleven-band tight-binding hamil¬ 
tonian (FTBH) obtained by Wannier transformation of 
density functional theory (DFT) results 28 . The trun¬ 
cated tight-binding hamiltonian (TBH) that retains only 
the first and partial second neighbor coupling is a good 
compromise, with little sacrifice of accuracy, while it 
yields clear interpretation of the physics in terms of the 
simple basis involved: (1) it contains all relevant orbitals 
near the Fermi level, which are p orbitals from X (chalco¬ 
gen) atoms and d orbitals from M (metal) atoms; (2) it 
preserves the phase and orbital information and character 
as well as all the crystal symmetries; (3) it preserves the 
orthogonality of basis functions which is crucial in con¬ 
structing interacting many-body theories; (4) it provides 
a systematic way of introducing higher order corrections 
and spin-orbit coupling terms; (5) it allows for interlayer 
hopping terms; (6) it allows the calculation of optical 
transitions and the Berry curvature of bands; (7) most 
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importantly, it contains no empirical parameters and is 
therefore a true ab-initio tight-binding hamiltonian. We 
also provide a simple recipe, based on correcting the band 
gap and band width, to show that scaling of the param¬ 
eters derived from DFT results gives a good representa¬ 
tion of the more accurate GW calculations. We only give 
these results for one case, MoS2, due to the considerable 
computational cost of obtaining GW band structure re¬ 
sults, but expect the procedure to be applicable to other 
similar materials. To augment the accuracy of the TBH, 
we also construct k • p hamiltonians for the low-energy 
bands around the T and K points of the Brillouin Zone 
(BZ) by expanding the FTBH. 

Building on the single-layer TBH and properly sym¬ 
metrized orbitals, we derive a set of hopping matrix ele¬ 
ments that accurately describe interlayer coupling within 
the tight-binding approximation. These matrix elements 
cover a range of possible distance between atomic sites 
in different layers, obtained by sliding two layers relative 
to each other. Surprisingly, the dependence of interlayer 
coupling parameters on distance shows very simple scal¬ 
ing with distance and is independent of angular depen¬ 
dence. This makes possible the calculation of the elec¬ 
tronic properties of configurations that involve arbitrary 
twists between successive layers, as well as calculation 
of inter-band optical transitions and Berry curvature of 
heterostructures consisting of different combinations of 
layers. Thus, our tight-binding scheme can be used to 
search efficiently for interesting behavior in a wide range 
of layer arrangements. 

The paper is organized as follows: In Sec. II we intro¬ 
duce the DFT and Wannier formalism as the tools em¬ 
ployed in studying the electronic structure of the relaxed 
TMDC crystal structures which are defined in Sec. III. 
To construct the TBH, we start with a monolayer unit 
in Sec. IV, where we discuss in detail the symmetries, 
spin-orbit coupling terms and the comparison to DFT re¬ 
sults. Based on the hamiltonian for the monolayer unit, 
we present in Sec. V the generalization to multiple lay¬ 
ers, modeled by introducing interlayer coupling through 
a transferable interaction term. In Sec. VI, we derive the 
complementary k • p hamiltonians, based on the relevant 
orbitals identified from the preceding analysis, that cap¬ 
ture all the low-energy physics of these materials. The 
applications for optical absorption and Berry phase with 
the TBH is discussed in Sec. VII. In Sec. VIII, we give 
a summary of the models and their potential applica¬ 
tions. The reduction of the TBH parameter space from 
symmetry considerations, the numerical values of the pa¬ 
rameters, and the GW calculations are presented in the 
Appendix. 


II. NUMERICAL METHODS: DFT AND 
WANNIER FORMALISM 

We perform the DFT calculations using the Vienna Ab 
initio Simulation Package (VASP) 29,30 . The interaction 


between ionic cores and valence electrons is described 
by pseudo-potentials of the Projector Augmented-Wave 
(PAW) type. The exchange-correlation energy of elec¬ 
trons is treated within the Generalized Gradient Ap¬ 
proximation (GGA) as parametrized by Perdew, Burke 
and Ernzerhof (PBE) 31 . We performed calculations with 
and without spin-orbit coupling to guide the construc¬ 
tion of the corresponding correction terms in the tight- 
binding hamiltonian. A slab geometry is employed to 
model single or double layers with a 20 A vacuum region 
between periodic images to minimize the interaction be¬ 
tween slabs. The crystal structure is relaxed until the 
Hellmann-Feynman forces are smaller in magnitude than 
O.OleV/A for each atom. The plane-wave energy cut¬ 
off we use is 450 eV with a reciprocal space grid of size 
25 x 25 x 1. Due to the underestimation of band-gaps 
and band widths in DFT, we also performed GW quasi¬ 
particle calculations to correct the M 0 S 2 DFT results. 

An alternative way to represent the DFT or GW band 
structure is to transform the Bloch basis into a basis 
of maximally-localized Wannier functions (MLWF) 32 as 
implemented in the Wannier90 code. To perform this 
Wannier unitary transformation for TMDCs, we use the 
seven highest valence bands and four lowest conduction 
bands composed of d orbitals for the metal atoms and 
p orbitals for the chalcogen atoms to form the localized 
Wannier functions. The initial projections are chosen to 
be the atomic p/d orbitals and the final converged Wan¬ 
nier functions are close to the localized atomic orbitals. 
The effective hamiltonian in the Wannier basis is inter¬ 
preted as the full-range ab-initio tight-binding hamilto¬ 
nian (FTBH). With a 25 x 25 x 1 fc-point grid sampling, 
the numerical accuracy of the FTBH is usually within a 
few meV compared to DFT or GW bands. Based on this 
FTBH, the truncated tight-binding hamiltonian (TBH) 
that retains only first and partial second neighbor cou¬ 
plings is often adequate to elucidate the nature of the 
bands. By adding more terms beyond nearest neighbor 
couplings, the numerical accuracy can be improved until 
it reaches DFT or GW-level results. This provides a sys¬ 
tematic way to improve the tight-binding hamiltonian. 

We wish to reiterate the importance of properly de¬ 
rived TBH parameters: the parameters obtained in the 
way discussed here are solely determined by the over¬ 
lap integrals for orbitals and the matrix elements from 
the full ab-initio calculation. The alternative approach 
in constructing a TBH often consists of optimizing the 
model parameters by fitting the band energies of ei¬ 
ther DFT calculations or experimental results. In that 
approach, the orbital character of the bands or the 
phases of the overlap integral are not explicitly consid¬ 
ered and maybe subject to overfitting. For example, in 
the nearest-neighbor tight-binding hamiltonian of single 
layer graphene, the band structure is invariant under a 
sign change of the nearest hopping parameter t. How¬ 
ever, the sign of this hopping parameter can be deter¬ 
mined from photoemission experiment 33 . We find that 
the sign we obtain from our procedure for constructing 
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FIG. 1. (a) Top view of a single layer TMDC crystal structure 
with hexagonal primitive vectors ai and a 2 . The smaller 
yellow circles represent the two chalcogen atoms separated 
by dx-x in the z direction and the larger brown circles the 
metal atoms, (b) Perspective side view for the TMDC 2H 
crystal structure composed of X A -M-X B monolayers with a 
the in-plane lattice constant and c/2 the separation between 
two units along the 5-axis. 

the TBH is consistent with Wannier analysis but is not 
determined from a blind fitting procedure. In some cases 
the discrepancy due to overfitting can be even more sub¬ 
tle. As another example, a study related to Bi-Sb alloy 34 
pointed out that it is possible to arrive at fitted hamil- 
tonians that are inconsistent in the description of sur¬ 
face states, which was attributed to the mistaken mirror 
Chern number. For topological materials, it is even more 
crucial to get the signs of the parameters right in order 
to give the correct topological invariants and Berry cur¬ 
vature. 


III. CRYSTAL STRUCTURE AND SYMMETRY 


sition and a three-fold rotation symmetry (AI 3 ); these 
symmetries are relevant to our hamiltonian construction. 
For multilayer TMDCs, there are various kinds of crystal 
structures, depending on the stacking along the c-axis, 
with periodic vector a 3 = cz where c is the distance be¬ 
tween repeating monolayer units. In this work, we focus 
on the 2H stacking with two monolayer units in a unit 
cell with hexagonal symmetry, shown in Fig. 1 (b), with 
the two units related by a mirror symmetry Ad 3 in the xz 
plane, or equivalently by a 7 r rotation around the 5-axis. 
In a monolayer unit or a 2H stack with odd number of 
layers, the inversion symmetry is explicitly broken while 
the symmetry is restored in a 2H stack of even number 
of layers. 

TABLE I. The in-plane lattice constant (a), unit cell size 
along the axis perpendicular to the plane (c), distance along 
the plane-normal direction between chalcogen layers (dx-x) 
and nearest neighbor bond between metal and chalcogen 
atoms (dx-M), all in A, for the TMDCs considered here. Num¬ 
bers in brackets and c are experimental bulk values 35 . 



MoS 2 

MoSe2 

ws 2 

WSe 2 

a(A) 

3.18 [3.16] 

3.32 [3.29] 

3.18 [3.15] 

3.32 [3.28] 

c(A) 

[12.29] 

[12.90] 

[12.32] 

[12.96] 

d X -x (A) 

3.13 [3.17] 

3.34 [3.33] 

3.14 [3.14] 

3.35 [3.34] 

dx-M(A) 

2.41 [2.42] 

2.54 [2.52] 

2.42 [2.40] 

2.55 [2.53] 


IV. MONOLAYER TIGHT-BINDING 
HAMILTONIAN 


TMDC crystals appear in different forms, labeled ac¬ 
cording to the stacking between successive layers as 1H, 
IT and IT’, the latter having a 2 x 1 reconstruction of 
the planar unit cell 10 . The electronic properties are tied 
to the underlying crystal structure. The semiconducting 
1H monolayer crystal structure is shown in Fig. 1 and 
is the focus of this work. The chalcogen atoms, labeled 
as X A /X B at the top/bottom, form a hexagonal lattice 
in each layer and project onto the same position in the 
plane of the middle layer of metal (M) atoms. As seen in 
Fig. 1(b), the M atoms have trigonal prismatic coordina¬ 
tion. The monolayer unit is characterized by the hexag¬ 
onal lattice constant a and the projected distance from 
X to the middle layer dx-x/2, which are given in Table 
I. In our convention, the primitive vectors are a\ = ax 
and a .2 = a(— \x + ^y) as shown in Fig. 1, with recip¬ 
rocal space vectors 61 = '^(x + ~^y) and b 2 = ~^^y- 
The BZ has special fc-points, F at the center, M= |foi, 
K ± =±1(263-62). 

The D 3 h symmetry group for the crystal contains a 
mirror symmetry in the xy plane (Adi), a mirror sym¬ 
metry in the yz plane (Ad 2 ) centered at each atomic po- 


A. Truncated Tight-Binding Hamiltonian (TBH) 

In the FTBH, the number of neighbor couplings in¬ 
cluded depends on the size of reciprocal space sampling 
in the Wannier construction. However, the couplings are 
short-ranged and fall off exponentially with the neighbor 
distance. Starting from a fully converged FTBH, we can 
truncate the hamiltonian to retain only the first-neighbor 
terms and partial second-neighbor terms to improve the 
accuracy. This TBH captures the essential features of 
the bands and their orbital character. To construct such 
a hamiltonian properly, we need to identify the relevant 
atomic orbitals and the symmetries of the crystal. For 
TMDC materials, the relevant atomic orbital basis for a 
monolayer are the d{p) orbitals of M(X) atoms, labeled 
as: 


Id = [4*.4y> d f x 2_ v 2 , dL, 4z. 




(1) 


The relevant symmetry operations of the monolayer in¬ 
clude the mirror symmetries Adi and Ad 2 and three-fold 
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TABLE II. Odd and even basis sectors (superscripts (o) and 
(e)) and their symmetries under xy mirror reflection (Mi) and 
yz mirror reflection (M 2 ) in terms of atomic Wannier orbitals: 
the d-like orbitals of M atoms and the p A /p B orbitals of X 
atoms located on the top/bottom layer. 


Index 

Basis Function 

Mi 

m 2 

1 


— d 

Uxz — ^xz 

- 

— 

2 

(°) 

Pz 

II 

s — y ^ 

^3 

- 

+ 

3 

= 71 (rf+rf) 

- 

+ 

4 

(o) 

Px 

= 71 (Px-Px) 

- 

- 

5 

( 0 ) 

Py 

= 75 ( P^-Py ) 

- 

+ 

6 


d^ e J = d z 2 

+ 

+ 

7 


II 

+ 

- 

8 

S 

X 

- d 2 2 

2_y2 — U x 2 -y 2 

+ 

+ 

9 

(e) 

Pz 

= 71 (rf-rf) 

+ 

+ 

10 

(e) 

Px 

= 7 ^(Px+Px) 

+ 

- 

11 

(e) 

Py 

= 7 ^Pt+Py) 

+ 

+ 


rotations 7 ^. 3 . They are used to classify the matrix ele¬ 
ments and reduce the number of independent parameters 
in the hamiltonian. First, the states are classified as odd 
or even under xy mirror symmetry M i 26 . The new basis 
that embodies this symmetry is: 



and contains five states in the odd sector and six in the 
even sector (see Table II for details). Under M 1 sym¬ 
metry the pair of X atoms in a unit cell is treated as a 
single composite atom with odd or even states. Because 
the hamiltonian commutes with Mi, it can be written 
in even and odd diagonal block form and there are no 
mixing terms between the two blocks which break the 
Mi symmetry. The two basis sets are linked by the uni¬ 
tary transformation T, where c/>| 0 = TV^. When there 
are multiple layers with interaction between layers, the 
transformation between the two sets of basis is needed 
to determine the inter-layer hopping terms. If spin is in¬ 
cluded, the generalized even/odd sectors of the Hilbert 
space can be defined (see the discussion on spin-orbit cou¬ 
pling). The minimal spinless hamiltonian for a monolayer 
is given by: 

^ (1L) = E^( fe ) JT S L) ( fc )^( fc ) ( 3 ) 

i,j,k 

This minimal tight-binding hamiltonian retains only 
the diagonal energy terms and the hopping terms as 
shown in Fig. 2: M-M hoppings (red) stand for the cou¬ 
pling from the odd (even) states on one M atom at to 



FIG. 2. Schematic diagram for all hopping terms in the 
TBH: M-M coupling (red) with six first-neighbor pairs; X-X 
coupling (green) with six first-neighbor pairs; X-M coupling 
(blue) with three first-neighbor pairs. Additional couplings 
between three second-neighbor X-M pairs (grey dashed) are 
used to improve the accuracy. The top bar indicates the la¬ 
bels of parameters that correspond to the different types of 
hopping terms, M' 1 , n = 1 ,..., 6. 

the odd (even) states on its six first-neighbor M atoms 
and X-X hoppings (green) are defined similarly for the six 
first-neighbor pairs of X atoms. X-M hoppings (blue) link 
states on M atoms to states at the three first-neighbor X 
atoms. The M-X hoppings are Hermitian conjugates of 
X-M hoppings. The grey dotted hoppings are additional 
X-M second neighbor terms that can be included to im¬ 
prove the accuracy for low-energy bands. All the rele¬ 
vant vectors that connect orbitals in the different hopping 
terms of the hamiltonian are defined in Table III. 


TABLE III. The hopping vectors in the TBH 



81 — ai, 82 — ai + a 2 , 83 — a 2 

#),#) 

d 4 = -(2ai + a 2 )/ 3, < 5 5 = (ai + 2a 2 )/3 
d 6 = (ai - a 2 )/ 3 

iW 

<^7 = —2(ai + 2a 2 )/3, 8 S = 2(2ai + a 2 )/3 
8 g = 2(a 2 - ai)/3 


We present next the matrix elements of the hamilto¬ 
nian after Fourier transformation to the k space in the 
BZ. We use the notation t\"j =< 4>i\H\(j)j > for the hop¬ 
ping matrix element from the state <pj to the state 4>i with 
the label s specifying the type and the distance between 
the two atoms where </>i and <ftj reside. The diagonal 
part for the hamiltonian contains e*, the on-site energy 
of orbitals, and the hopping terms between orbitals of 
the same type at first-neighbor positions. The diagonal 
terms of the tight-binding hamiltonian take the form: 

H\f\k) = ei + 2tf}cosk-8i 

+ 2 tfj [cos(fe- 8 2 ) + cos(fe- <S 3 )] 


( 4 ) 
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The off-diagonal hopping matrix elements between the 
same type of atoms (M-M and X-X) can be classified 
into two categories depending on the symmetry of the i 
and j orbitals under M 2 ' for (*,J)=(3,5),(6,8),(9,11), the 
symmetry is (+), giving 

(fc) = 2 tg cos k- 8\ + tg [e~ ik S > + e^] 

+ tf][e ikS2 +e ik ' S3 ] 

while for (*j)=(l,2),(3,4),(4,5),(6,7),(7,8),(9,10),(10,11), 
the symmetry is (—), giving 


H[f\k) = - 2iig sin k- 8 1 + tf][e~ ikS2 - e ~ ikS3 ] 


+ t\ 3 -[-e lkS2 


i,3 

Jk-5 3 ] 


( 6 ) 


A different type of hopping connects M and X atoms. 
Each M atom has three first neighbor X pairs. The ma¬ 
trix elements can still be classified by the symmetry of or¬ 
bitals i and j under M 2 , that is, for the pairs (*j)=(3,l), 
(5,1), (4,2), (10,6), (9,7), (11,7), (10,8), the symmetry is 
(+), giving 


Hlf\k)=tf ) j [e iks *-e ik - s °] (7) 

while for the pairs («,))=(4,1), (3,2), (5,2), (9,6), (11,6), 
(10,7), (9,8), (11,8), the symmetry is (-), giving 

tfg L) ( k ) = tg [e ikSi + e ikS& ] + tf]e ikS5 (8) 

The hermitian character of the hamiltonian requires 
H^\k) = Hj lk \k)* and otherwise unassigned 

H\^f\k) terms are zero. 

The above minimal tight-binding hamiltonian has 86 
parameters for all ej, tg, tg, tg(, tg and tg under 
the Mi and .M 2 symmetry classifications. The three¬ 
fold symmetry 77-3 can be used to further reduce the size 
of the parameter space by identifying equivalent coupling 
directions: with this simplification, only a subset of 
and all tg, tg are independent parameters while tg, 

tg and ig can be expressed as linear combinations of 
those, which reduces the number of free parameters to 36. 
The relationships that express these symmetry-imposed 
simplifications are given in Appendix A. 

In most applications, the only relevant degrees of free¬ 
dom involved in physical processes are the lowest con¬ 
duction band and highest valence band at K. The orbital 
character of these bands are given in Table IV. They are 
even states under Mi. 

The bands in the above minimal tight-binding hamil¬ 
tonian can be improved in an efficient way by includ¬ 
ing terms beyond first-neighbor couplings. Based on an 
analysis of the orbital character of low-energy bands (see 
next section), our choice is to add the dominant second- 
neighbor X-M coupling terms between the even states 




M K 


r r m k 


FIG. 3. TBH (red lines) band structure along T-M-K-r and 
the comparison with DFT results (blue circles) for monolayer 
spinless (a) M 0 S 2 (b) WSe 2 - 


(described by the grey dashed lines in Fig. 2). There are 
four additional independent parameters in H ^ 1L ) needed 
to include these coupling terms, labeled tg, 4i4> 4% 
fg g , and the correction to the TBH reads: 

H'^ = ^\{k)H\f\k)^{k) (9) 

i,j,k 

with the hermitian matrix elements: 


J?gS L) (fc) = t ( 9 6 l(e lkS7 + e ikSa + e lkS9 ) 
H[[f(k) = tg 6 (e ifc ^ - ±e ik '- \e ikS *) 

H^{k) = ^fg 6 (-e ifc -*> +e ikS ») 

H a f\k) = £% ik * - \e ik - s * - l -e lk5 *) 

H'<>f\k)='ft ( gl(-e ikSa +e ikS °) 

H'$\k)=^t%(e ikSa +e ikS °) 

H$\k) = H$\k) = ^tfl 8 {e ik Sa - e ik S9 ) 

= 4?8(e ifc ' aT + + \e lkS °) 


The final TBH is defined as the sum of the two terms, 
^(il.tbh) = #(iL) + 7 j'(il)_ In Fig . we show the 

resulting band structure from the TBH for M 0 S 2 and 
WSe 2 in comparison with DFT results. We also note 
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that even when the hamiltonian is truncated to retain 
only the first-neighbor terms (only H U L )), this minimal 
hamiltonian still preserves the essential features and or¬ 
bital character of the bands. The tight-binding results 
for the valence bands are better than for the conduction 
bands. This is expected due to the proximity of other 
higher energy bands near the conduction bands to which 
we have restricted the model presented here. If the bands 
belong to the same symmetry group representation, the 
level anticrossing between these bands introduces orbital 
character mixing. This mixing generates higher order 
corrections in the effective hamiltonian when the irrel¬ 
evant bands are integrated out. To account for these 
effects, higher order coupling terms which involve cou¬ 
pling between farther atoms are needed to obtain a more 
accurate description. 

To improve the underestimated band-gap value at the 
DFT level, we performed GW calculations and corrected 
the tight-binding hamiltonian parameters accordingly for 
M 0 S 2 . The details and the scaling factors of the tight- 
binding parameters for the GW quasi-particle energies 


are given in Appendix B. For the rest of the discussion, 
we will focus on the TBH bands based on DFT results, 
which is adequate to bring out the physics. 

B. Orbital Characters of Bands 

It is interesting to investigate the spinless orbital 
character of the bands at the high-symmetry F and 
K± points, which are relevant for low-energy dynamics. 
These k points have 77-3 symmetry, and the bands can 
be characterized by the symmetry group representations. 
To facilitate the classification of the states, we assign as 
labels the eigenvalue m z of the angular momentum opera- 

tor L z . d~ j-2 ( d x 2 _y2 ziz id X y ), d ~si -F ( d xz A idy Z ), 

do = d z 2 , p ±1 = T^(Px ± Wy), Po = Pz- Under clock¬ 
wise 77 3 symmetry with the axis of rotation centered at 
the atomic positions, these orbitals are invariant and have 
eigenvalues e + z27rm U3 i The orbital character of the bands 
at these special k points is given in Table IV, along with 
the coefficients Ci ( Ci = \J\ — cf) for the four TMDC 
materials considered here. 


TABLE IV. Orbital composition and numerical coefficients for the eleven bands at T, K+ (states at K_ are related to the states 
at K+ by the time-reversal symmetry). In the subscript notation fh stands for — m (m =1,2). For the coefficients, we define 
Ci = 1 — c?. Bands 1-7 (8-11) are valence (conduction) bands. 


K+ 


M0S2 

MoSe 2 

ic 8 d[ 0> - c 8 p l 0 o> 

cs 

0.5263 

0.5772 

c 4 dj 0> + c 4 p[ o) 

c 7 

0.4255 

0.4344 

c 7 <% ] - c 7 p^ e) 

Co 

0.4432 

0.4204 

ic 5 d ( 0 c) + c 5 pj C> 

C 5 

0.4026 

0.4012 

ic 6 d. 2 c> + c 6 Pi 

c 4 

0.6268 

0.6149 

c 8 d[ 0> - icsp^ 


ws 2 

WSe 2 

[°1 

P) 

cs 

0.4826 

0.5394 

C 7 df + c 7 p [ 0 c) 

c 7 

0.3883 

0.3988 

c 6 d ( 2 } + ic 6 p[ c> 

Co 

0.4643 

0.4450 

c 5 dl Q) + ic 5 p\ e) 

^5 

0.3564 

0.3536 

c 4 d ( ° ] — c 4 p{ } 

c 4 

0.6291 

0.6173 


Band 

r 


11 

ic 2 d 2 c> - c 2 p\ c> 

10 

ic 2 d t eJ - c 2 Pi 

9 

c 3 c4 0) + C 3 p[ 0> 

8 

c 3 d\ 0> + c 3 p\° } 


MoS 2 

MoSe 2 

7 

Ci dp - ciPq 

c 3 

0.7239 

0.7779 

6 

c 3 d[ 0> - c 3 p[ 0> 

c 2 

0.8010 

0.8234 

5 

c 3 d\° ] - c 3 p\ o) 

Cl 

0.5711 

0.5348 

4 

(°j 

Po 


WS 2 

WSe 2 

3 

c 2 d 2 c) - ic 2 p\ c> 

C 3 

0.6698 

0.7322 

2 

c 2 d\ ) e) - ic 2 pi 

c 2 

0.7850 

0.8077 

1 

cid\ c> + cip^ } 

Cl 

0.5729 

0.5414 


At T, there are four pairs of doubly degenerate states 
which form the two-dimensional representation under the 
symmetry group of the crystal. Each of the states are 
chosen to be invariant under 1Z 3 and the symmetry de¬ 
mands m p — md = 3 n for the state mixing, where m is 
the eigenvalue for the orbital angular momentum L z and 
n is an integer. Under yz mirror symmetry, m goes to 
—rn, and this results in the two dimensional degenerate 
subspace at T. 

At K-t, all states are non-degenerate. Because the 
wavevector k at K-t is not zero, we have to consider the 


symmetry transformations of the orbitals and of the spa¬ 
tial wavefunction. We take the rotational axis of 1Z 3 to 
be located at the position of M atoms. Since X atoms are 
not at the center, they will acquire an additional phase 
factor, e ± * 27r / 3 , under this symmetry. When both phase 
factors from orbital and spatial transformations are in¬ 
cluded, we arrive at the constraint m v ± 1 — ma = 3 n! for 
each state and the eigenvalues of 'R. 3 can be determined. 
States at K_ are related to states at K + by complex con¬ 
jugation of the wavefunctions in the ij) basis, dictated by 
time-reversal symmetry for spinless particles. 
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As an example of the implications of these symmetry 
considerations, we mention that the highest valence band 
and lowest conduction band at K-t are coupled by pho¬ 
tons of polarization a± due to the broken inversion sym¬ 
metry and the chiral selection rule (see the application to 
optical transitions in Sec. VII). This allows the optical 
control of valley degrees of freedom and the breaking of 
valley degeneracy by circularly polarized light 14 . The it 
light with polarization perpendicular to the layer, which 
has odd symmetry under Mi, cannot couple these two 
even states 36 . In forming stacks of TMDC layers, the 
bands with strong po character should be affected most 
by the interlayer hybridization. This implies that the 
highest valence band at T is sensitive to the interlayer 
coupling. On the other hand, the low-energy bands at 
K± are not sensitive to interlayer couplings because they 
are dominated by the d orbitals of M atoms, which are 
deeper inside the monolayer unit. 


C. Spin-Orbit coupling 

Up to this point, we have not included the spin degrees 
of freedom. Due to the lack of inversion symmetry in 
TMDCs, spin-splitting of bands generally occurs over the 
entire BZ. However, the bands are still doubly degenerate 
at the time-reversal invariant points T and M due to the 
Kramers’ theorem and at special k points, such as along 
the T-M direction, which are protected by crystal mir¬ 
ror symmetry. In TMDCs, the low-energy bands at K± 
are not spin-degenerate and hence the spin is crucial for 
understanding the low-energy dynamics. The even/odd 
sectors under Mi mirror symmetry must be redefined to 
incorporate the spin degrees of freedom. Under Mi xy 
mirror operation, S will flip the signs of S x and S y com¬ 
ponents but not of S z . Hence, Mi mirror symmetry acts 
as a 7r rotation along the 0 -axis, Om, = e l ' KSz ^ h . The 
spin up (down) state is an eigenstates with eigenvalue 34 
+i (—i). The generalized Mi operation is the product 
of its action on the orbital and on the spin degrees of 
freedom. In this extended Hilbert space, the states of 
even orbitals with spin up (down) and the states of odd 
orbitals with spin down (up) are in the same enlarged +* 
(— i) sector under M\. These generalized ±i sectors dic¬ 
tate which set of states can be mixed in the hamiltonian 
and there are no mixing terms between the generalized 
+i and —i sectors. 

To incorporate this spin-orbit coupling effect in the 
tight-binding hamiltonian, we approximate the contri¬ 
bution for spin-orbit interaction by the atomic term 
X^q X L-S for each individual M and X atom 3 '. In this 
expression, L = L x , L yi L Z (S = S x , S y , S z ) stands for the 
orbital (spin) angular momentum operators, with A so 
the strength for the spin-orbit coupling term which de¬ 
pends on the atomic species (see Table VIII). We obtain 
the values for various M and X atoms from the energy 
splitting of a single atom in a large unit cell as calcu¬ 



FIG. 4. Spin splitting for the highest pair of valence bands 
in M 0 S 2 : (a) along the T-K direction as obtained from TBH 
+ SOC (red lines) and DFT (blue dots) calculations (b) The 
angular dependence with fixed k\\=K/2 and fcy=K of TBH + 
SOC (red dots); the black lines are the fitted cos(3$) func¬ 
tions. 


lated by DFT when spin-orbit coupling is included in the 
Kohn-Sham equations. The Hilbert space for the mono- 
layer hamiltonian with spin-orbit coupling (SOC) is the 
direct product space of Wannier orbitals and spin de¬ 
grees of freedom, with 22 states in total. The TBH + 
SOC hamiltonian is given by: 


#so L) = E[ 4 ( fe ) jff tt L ’ TBH) ( fc )^t(fc) 

( 11 ) 

+ ^{(fc)^ L ’ TBH) (fc)^(fc) + ft(k)H hS $(k)\ 

The diagonal blocks in the first term Iuj;j L ' TBH ' = 

= ij( 1L TBH) are tbjj described above. 
These are the spin-independent hopping processes. The 
effect of spin-orbit coupling, Hls, is incorporated by the 
on-site A soL- S term for each atom. Because it is an on¬ 
site term, it does not carry momentum dependence and 
is a constant matrix with the matrix elements: 


=<4,a|(A^oi M + Af 0 ^ + Af 0 L|)-5|^> 

It is straightforward to evaluate these matrix elements 
with the help of Clebsch-Gordan coefficients. The LS 
terms can also be written as ( L + S _ + L~S + )/2 + L Z S Z , 
and the change in the evenness/oddness of the orbitals 
is always accompanied by the change of the spin state. 
This observation affirms the above classification of the 
generalized +i and — i sectors under Mi. 

To compare the TBH + SOC results with the full 
DFT results, we compute the spin splitting energy, Aj£ = 
e|(fc) — ej (fc), for the two highest valence bands which 
are a spin up and down pair. The spins remain dou¬ 
bly degenerate along the T-M direction and the largest 
splitting is along the T-K direction. In Fig. 4 (a) the 
TBH + SOC spin splitting is plotted along T-K (in red 
line) and compared with the DFT results (in blue dots). 





FIG. 5. V PPl o(r) (upper) and V pp ^(r) (lower) for interlayer 
pairs as function of the pair distance r in a M 0 S 2 bilayer unit. 
Black, red, green, blue, cyan, brown are data points from 1st, 
2nd, 3rd, 4th, 5th, 6th nearest distances. The solid black lines 
going through the points are the fits as obtained from Eq. 17. 
The inset is the perspective view for the interlayer coupling 
between one of the S atom on the top unit (orange) to all 
other S atoms on the bottom unit (yellow). 


The change in the spin splitting is also accompanied by 
the change in the underlying orbital composition from 
T to K 38,39 . Because of the crystal symmetry, this spin 
splitting is expected to have the form cos(30) with 
out-of-plane polarization 38 . In Fig. 4 (b), we compare 
the angular dependence of the spin splitting at fey = K/2 
and fcy = K in TBH + SOC (red dots) with the fitted 
cos(30) formula (black line). This fitted formula holds 
well over the BZ and the spin splitting for the BZ can be 
determined with the radial data in Fig. 4 (a) multiplied 
by the angular dependence, that is, the cos(3 9) term. 
The radial part at small k rises as k 3 , consistent with the 
dominant spin splitting term ( k x + ik y ) 3 + ( k x — ik y ) 3 
from the crystal symmetry. A more accurate SOC effect 
will be presented in connection to the k • p hamiltoni- 
ans at the K-t points which are tailored for low-energy 
physics with higher accuracy. 


V. INTERLAYER PAIR COUPLING 


present in the exchange-correlation functional, but is not 
taken explicitly into consideration here. The van der 
Waals force affects mostly the crystal structure and the 
equilibrium vertical separation of the layers. The 2H 
crystal structure is shown in Fig. 1 (b) and the equi¬ 
librium distance of layers along the c-axis is taken from 
experimental values. 

The dominant interaction contribution to the elec¬ 
tronic band structure comes from orbital hybridization 
between adjacent X layers. These are the p-p hop¬ 
ping terms between X atoms, with p z -p z coupling be¬ 
ing the dominant contribution due to the orbital orien¬ 
tation. We will show that the two-center Slater-Koster 
approximation 40 describes interlayer coupling well and 
transferable empirical interactions can be obtained for 
the strength of a and 7 r bonds as functions of the pairing 
distance across the spacing between layers. 

To extract the interlayer interaction, the full-range 
Wannier tight-binding hamiltonian is constructed for a 
2H bilayer TMDC unit as in the monolayer case 41 . The 
initial orbital projections are the same atomic orbitals 
and the Hilbert space twice as large. The hamiltonian 
can be separated into two parts: the first corresponds to 
the coupling within the same monolayer unit and it is in 
diagonal blocks of each layer, while the second, H> nt , 
includes the interlayer hopping matrix elements, fy b) : 

^int L) = H ( r2 )4X ( r2 “ ri )^!.P J ( ri ) + hx - 

P'n r 2,Pj,ri 

(13) 

For each interlayer pair of X atoms, nine hopping pa¬ 
rameters are needed to describe the interactions between 
three p orbitals on each, denoted as: 

t ( p ^l O 2 - ri) =< 0 2 , P ' (r 2 )|4nt L) l^i , Pi ( r i) > ( 14 ) 

Within the two-center Slater-Koster approximation 40 , 
these nine p-p coupling parameters are reducible and can 
be simplified to two independent parameters, V PPjCr and 
Vpp,Tn corresponding the a and 7r bonds. The nine hop¬ 
ping terms for each pair reconstructed from the V PPtCr and 
Vpp^ representation are given by: 


A. p-p Interlayer Coupling 

In multilayer structures or the heterostructures of sev¬ 
eral different 2D layers, attractive forces of van der Waals 
type bind the individual layers together. The resulting 
electronic states show interesting new features due to the 
interlayer interactions. In this section we investigate the 
induced changes in electronic band structure and the in¬ 
teraction hamiltonian between layers of the same kind. 
We stress that the many-body physics of the van der 
Waals interactions is not included in the electronic band 
structure. In principle, this many-body effect should be 


#2 M = ( v rpAr) ~ V PP Ar)) ^ + Vp P Ar)ki (15) 

with r = |r|. The inverse relations are used to extract 
Vpp,<r(r) and V pp ^{r) from these nine hopping parameters 
for each interlayer pair: 



Vpp,a(r) — 


(r) 


nr 3 

J.2 


( 16 ) 
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For a bilayer configuration, multiple pairs of interlayer 
distances can be identified starting from the dominant 
contribution of the shortest-distance pair. These are in¬ 
dicated by the arrows in different colors in the inset of 
Fig. 5. 


B. Universal Form for Interlayer Coupling 

A single fixed configuration of 2H bilayer TMDC unit 
gives pairs with distances at only a discrete set of values. 
A spatial translation or a twist in the crystal orientation 
for each layer unit can exist in the real heterostructure. 
To study the interlayer coupling in a more general crys¬ 
tal structure, we need to map out V pp i(T and V pp j7r as 
functions of the pair distance r treated as a continuous 
variable by varying the conventional 2H bilayer structure. 
We use a horizontal translation of the top unit relative to 
the bottom unit of the conventional 2H bilayer to obtain 
interactions as the distance of the pair is varied contin¬ 
uously. The vertical distance between the layers, c/2, is 
kept fixed when the relative translation is applied. Since 
there is no rotation between the two monolayer units, 
the primitive cell is the same with displaced atoms at 
the top layer. The vector for the translation is spanned 
by l\CLi + 12&2 with li between 0 and 1. For each crys¬ 
tal configuration that includes a translation we extract 
the sets of interlayer pairs and the V pp ^ a and V pp ^ values 
associated with them. As shown in Fig. 5, V PPtCr and 
Vpp,ir values from the interlayer pairs in various crystal 
configurations collapse on two single curves respectively. 
This demonstrates the validity of the two-center Slater- 
Koster approximation and the pair distance r as the only 
variable needed to parametrize the interlayer pair interac¬ 
tion. We use an exponential form to fit V PP:b as functions 
of r where b=a, 7 r . 42 

Vp P ,b{r) = v b e*p(-(r/R b ) rib ) (17) 

This interlayer interaction depends only on the species 
of X atoms, since M atoms are well inside the mono- 
layer unit and their orbitals do not contribute directly to 
the interlayer interactions. We give the values of these 
parameters in Table V for S-S and Se-Se interlayer inter¬ 
actions. 


TABLE V. Interlayer V vv , a and V vv ,-n parameters. 


Parameters 

V s-S 

PP,& 

yS-S 

r PP,7T 

T/-Se—Se 
V pp,cT 

ySe—Se 
v pp,iv 

v (eV) 

2.627 

-0.708 

2.559 

-1.006 

R( A) 

3.128 

2.923 

3.337 

2.927 

V 

3.859 

5.724 

4.114 

5.185 


The validity of the two-center approximation can be 
attributed to the weak coupling of the two units in the 
bilayer structure. The Wannier functions of p orbitals at 


X atoms show only small crystal field distortion. That 
is, the atomic X orbitals are insensitive to the local crys¬ 
tal environment, or the crystal orientation. The inter¬ 
layer coupling parameters do not require an angular de¬ 
pendence as established by the results shown in Fig. 5. 
With this azimuthal symmetry, the pair distance infor¬ 
mation is enough to determine the interaction strength, 
irrespective of the pair orientation. This also means that 
the interlayer interaction will have the same form when 
the two units are rotated by an arbitrary angle relative 
to each other, which makes our TBH scheme applicable 
to heterostructure that involve a twist between layers. 

The azimuthal symmetry of atomic orbitals that par¬ 
ticipate in the interlayer coupling is not generally ex¬ 
pected due to the crystal field distortion of atomic or¬ 
bitals from the presence of neighboring orbitals. One 
salient case is the tight-binding hamiltonian for bilayer 
graphene 43 in Bernal stacking. Two types of pairs, 73 
and 74, share the same pair distance but have very dif¬ 
ferent hopping strength. In general, the crystal field in¬ 
troduces higher angular momentum mixing to atomic or¬ 
bitals. This gives rise to additional orientational depen¬ 
dence for interlayer interactions beyond the simple inter¬ 
layer pair distance. A more general interaction model 
can be constructed to account for both pair distance and 
orientation, which are crucial for modeling generic two- 
dimensional layered materials. 


C. Application to 2 H Bilayer and Bulk 

To test the interlayer interaction hamiltonian and com¬ 
pare with the full DFT calculation, we construct the 
tight-binding hamiltonian in Fourier space for a M 0 S 2 
2H bilayer as follows: 

H (2h) = £[# (k)H[ lh \k)Mk) + fyk)H' 2 {1L \k)Mk) 

k 

+ $l(h)v£ L \k)Mk) + h.c.} 

(18) 

Within each monolayer unit, i/( TBH ) is used for the 
hamiltonian, H^ 1L \ The relative rotation for the sec¬ 
ond monolayer unit, H 2 1L \ can be taken into account 
by either a rotation in the hamiltonian basis or xz mirror 
operation 23 . For the interlayer hopping V^ L \ all inter¬ 
layer pairs with distance less than 5A are included. Due 
to the extension of d z 2 orbitals into the interlayer region, 
additional interlayer coupling between d z 2 and p z can be 
included to improve the accuracy of the highest valence 
bands at T which are composed of d z 2 and P 2 orbitals. 
These interlayer couplings are 60 meV for the nearest X- 
M and 26 meV for the second nearest X-M pairs. We 
show in Fig. 6 (a) the band structure with this bilayer 
tight-binding hamiltonian (red lines), which is in good 
agreement with the full DFT results (blue circles). 
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FIG. 6. (a)TBH band structure (red lines) along T-M-K- 
T and the comparison with DFT results (blue circles) for a 
M 0 S 2 bilayer in the 2H structure, (b) The overlap matrix ele¬ 
ment strength between two layers with the interlayer coupling 
hamiltonian that leads to splitting of single-layer states, on a 
scale from red (strongest) to white (weakest). 


To gain insights of the effect of interlayer interactions, 
we can treat the interlayer couplings W l l) in the context 
of a perturbative hamiltonian. When this perturbation 
is set to zero, the bands for the bilayer unit are dou¬ 
bly degenerate. For a generic k point, we choose the 
eigenstates, ipi (k) and <p 2 {k) to exist on the individual 
layers respectively. The degenerate subspace is spanned 
by these two states. When the interlayer coupling V ^ LL - ) 
is turned on, it mixes the two states by introducing off- 
diagonal matrix elements, < </Ji(fc)|W LL )(fc)|</? 2 (fc) >• In 
Fig. 6(b), the M 0 S 2 monolayer bands are plotted on a 
color scale that represents the magnitude of the inter¬ 
layer off-diagonal matrix element with red representing 
the largest values and white the smallest values. The 
splitting between two degenerate states is proportional 
to this overlap matrix element. The strongest coupling 
and splitting comes from the bands with dominant p z 
character near the T point. This hybridization drives the 
direct band gap into an indirect one when more than one 
layer is included in the structure. 

The interaction does not always guarantee the off- 
diagonal mixing between these two states due to sym¬ 
metry selection rules. The symmetry and the orbital 
character (see Table IV) can be used together to explain 
the degeneracies at K for the 2H bilayer. The degener¬ 
acy occurs when the two states, <^>1 (K) and </? 2 (K)i at the 
same energy transform differently under IZ 3 symmetry. 
As a result, the off-diagonal matrix elements of inter¬ 
layer coupling vanish in this case. Due to the restoration 
of inversion symmetry in the 2H bilayer, the bands are 
doubly degenerate with two spins when spin-orbit cou¬ 



r m k rr' M' k 1 r 


FIG. 7. Band structure of the bulk, 3D crystal of M 0 S 2 , 
in the 2H structure, as obtained from DFT calculations (blue 
circles) and from the TBH (red lines) derived from the bilayer 
system (including inter-layer coupling). The symbols T-M-K- 
T correspond to k z = 0, while those with prime correspond 
to k, = it/c. 

pling is included. 

For the bulk band structure in the 2H configuration, 
there will be additional interlayer hopping terms under 
periodic boundary conditions compared to the bilayer 
unit. To have the minimal change to the monolayer 
hamiltonian and to incorporate properly the k z depen¬ 
dence, we make the gauge choice that for the 2 direction 
atoms in the same layer are treated as if they are on the 
same z plane. All the k z dependence enters into the inter¬ 
layer terms and effectively involves a shift by c/2 in the z 
direction for the interlayer pair. This gauge choice makes 
it easier to write the hamiltonian within each monolayer 
unit and renders them the same as the case with k z = 0. 
The resulting bulk bands from the tight-binding hamilto¬ 
nian (red lines) are shown in Fig. 7; these results clearly 
show that the bulk band structure is in good agreement 
with the full DFT results (blue dots). 

VI. K P HAMILTONIAN AT BAND EXTREMA 

For the low-energy excitations around the band gap 
edges, it is useful to have an accurate perturbative ap¬ 
proach in the form of a k • p hamiltonian. We construct 
this type of hamiltonian by performing a series expan¬ 
sion of the FTBH. In contrast to the TBH which faith¬ 
fully represents the symmetry, band structure and orbital 
character over the entire BZ, the k • p hamiltonian is con¬ 
structed and optimized for small regions of the BZ with 
selected bands where the low-energy excitations reside. 
We perform this k • p expansion for: (1) the highest 
valence band at T; (2) the lowest conduction and the 
highest valence bands at K±. The expansions are based 
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on the hamiltonian of the spinless electrons of the sin¬ 
gle layer for each material. The corresponding correction 
terms for the spin-orbit coupling or the interlayer inter¬ 
action in the bilayer structure can be introduced as ad¬ 
ditional terms. In general, the presence of these pertur¬ 
bation terms would renormalize the parameters for the 
unperturbed hamiltonian of the single layer. We include 
these higher order effects as needed, and keep only the 
necessary lowest order terms otherwise. 

The derivation of the k • p hamiltonian from the se¬ 
ries expansion of the FTBH, as opposed to fitting the 
band structure of ab-initio results 13 ’ 25 , allows us to keep 
the phase and the orbital information from the analysis of 
DFT and FTBH results. The spinless k • p hamiltonian of 
the single layer is expanded up to second order in fc-space 
to capture the effective masses for the bands. There are 
two contributions for the second order terms, originat¬ 
ing from the second order expansion of the hamiltonian 
within the subspace or the second order virtual transi¬ 
tion terms. The latter comes from the effective terms by 
integrating out the irrelevant bands through a Schrieffer- 
Wolff transformation 44 . 

The basis refers to the Bloch bands at T and K± and 
consists of hybrids of p and d orbitals as given in Table 
IV. At the r point, the valence band is of predominantly 
d z 2 character. For the K-|- point, the conduction band is 
mostly of d z 2 character while the valence band is mainly 
composed of ( d x 2_ y 2 + ird xy ) / y /2 orbitals with r = ±1 
for K±. The gauge choice for the basis is to have as the 
real positive component either d z 2 or d x 2 _ y 2 , depending 
on which one dominates. The group of the wave vec¬ 
tor dictates the form for the hamiltonian under symme¬ 
try transformations 45 . Time reversal symmetry links the 
hamiltonian at K±. In the following discussion, er acts 
on the conduction and valence degrees of freedom at K 
and a x , a y , a z are the Pauli matrices, s acts on the spin 
subspace while /t acts on the layer index when a bilayer is 
considered. The direct product of these operators and the 
Hilbert space is implicitly assumed when extra degrees of 
freedom and the perturbative terms are included. 

We begin with the k • p spinless hamiltonian for a single 
layer. The zero energy point is the highest valence band 
at K and k = ( k x ,k y ) is the relative crystal momentum 
from the expansion point. At the F point, the valence 
band is an isotropic parabolic band due to IZ 3 symmetry. 
The parameter g 0 is the relative energy shift of valence 
bands with respect to K and a the lattice constant. These 
considerations lead to the following form: 


di F (k) = g 0 +gia 2 \k \' 2 (19) 

The structure of the k • p hamiltonian at K-t is 
richer because both conduction and valence bands are in¬ 
volved. Time-reversal symmetry demands H- T (k x , k y ) = 
H*(—k x ,—k y ) between two valleys and the hamiltonian 
reads 24 



FIG. 8. (a) Comparison of the k • p hamiltonian (red lines) 
and the DFT results (blue dots) for the M 0 S 2 monolayer with 
spin-orbit coupling at K+. (b) Energy contours in eV for the 
highest valence band centered at K+. The deviation from 
spherical symmetry is due to the term. 


y.^(k) = y (1 + <J Z ) + f\a{rk x G x + k y d y ) 

+ a 2 |fc| 2 (/ 2 + f 3 a z ) + f±a 2 (k 2 x d x - k 2 a x - 2 rk x k y a y ) 

( 20 ) 

In general, other off-diagonal forms which are consistent 
with the crystal symmetry can also appear in a different 
gauge choice of the basis. Up to linear order, this hamil¬ 
tonian has the form of the massive Dirac equation with 
energy gap /o. Expanding this 2x2 hamiltonian, the 
eigenvalues and the dispersion up to second order for the 
conduction and valence bands are approximately given 
by: 


e c (k) ~ /o + (/*2 + h + ^) a2 |&l“ 
e u (fc)«(/2-/ 3 -^)a 2 |fc| 2 


( 21 ) 


The effective masses are obtained from the coefficients of 
the quadratic terms and are dominated by the massive 
Dirac terms / 2 // 0 while the / 3 (/ 2 ) term gives additional 
particle-hole symmetric (asymmetric) effective mass con¬ 
tributions. 

When spin is included, the Hilbert space is doubled 
and there will be additional spin-orbit coupling terms. 
At the T point, the bands are doubly degenerate due 
to Kramers’ theorem and there are no additional spin- 
orbit coupling terms to zeroth order. The spin splitting 
is non-zero when moving away from T and it scales as fc 3 
(see the discussion for Fig. 4). At the K point, spin-orbit 
coupling splits the two spin states for both the conduction 
and valence bands. The lowest order four-dimensional 
spin-orbit coupling term is: 


A H* 


hr{ 


1 - <7, 


)s z + ff>r( 


1 + d z 


)Sz 


( 22 ) 


2 


2 
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where (fo) represents the spin-orbit coupling strength 
for the valence (conduction) band. The small spin¬ 
splitting term / 6 for the conduction bands at K is more 
subtle from the two contributions in the perturbative 
expansion 4 ' 5 . The first part is the negative contribution 
from the projection of the spin-orbit coupling operator 
within the conduction band subspace, while the second 
positive contribution is from the virtual coupling to other 
higher bands by the spin-orbit coupling. These two con¬ 
tributions tend to cancel out and result in the alternat¬ 
ing signs in f e in MoX 2 and WX 2 as reported in the 
literature 4 ' 1 . Due to time-reversal symmetry and broken 
inversion symmetry, the spin state flips when changing 
valley and the combination ts z can be viewed as the ef¬ 
fective coupling term between valley and spin. In the 
presence of spin splitting, the gap between bands is dif¬ 
ferent and this will affect the dispersion compared to the 
spinless case. The dispersion for the bands with spin is 
given by the same expression as in Eq. (21) with the 
gap fo replaced by the new spin-valley dependent gap 
= fo+Ts z (f 6 - f 5 ). 

In Fig. 8 (a), the MoS 2 k ■ p hamiltonian is compared 
to the DFT bands around the K + point. In Fig. 8 (b), we 
plot the energy contours for the valence bands at several 
energies by computing the spectral function A(k,ui) = 
—Im[Tr[(w — 1-L + *<5 ) -1 ]]/n with a smearing 5 = 2 meV. 
The / 4 term gives the deviation of the energy contours 
from spherical symmetry due to the underlying three-fold 
crystal symmetry 4 '. 


TABLE VI. k ■ p hamiltonian parameters in eV with the lat¬ 
tice constant a in A. 



MoS 2 

MoS 2 GW 

MoSe 2 

ws 2 

WSe 2 

a 

3.18 

3.18 

3.32 

3.18 

3.32 

go 

—0.0167 

-0.1161 

-0.2712 

-0.0648 

-0.3347 

91 

-0.1173 

-0.1004 

-0.0642 

-0.1314 

-0.0680 

92 

0.1563 

- 

-0.0705 

0.0940 

-0.1570 

93 

-0.2106 

- 

-0.1351 

-0.2365 

-0.1534 

9i 

0.1699 

- 

0.1464 

0.1869 

0.1599 

95 

-0.3319 

- 

-0.3352 

-0.3272 

-0.3205 

fo 

1.6735 

2.4826 

1.4415 

1.8126 

1.5455 

h 

1.1518 

1.7381 

0.9560 

1.4073 

1.1894 

/2 

0.0744 

0.0957 

0.0494 

0.1551 

0.1184 

fs 

-0.0613 

-0.2917 

-0.0493 

-0.0175 

-0.0064 

h 

-0.0780 

-0.1308 

-0.0654 

-0.0709 

-0.0627 

/*5 

0.0746 

- 

0.0929 

0.2153 

0.2335 

h 

-0.0015 

- 

-0.0106 

0.0148 

0.0180 

h 

-0.0417 

- 

-0.0508 

-0.0502 

-0.0612 


In the 2H bilayer unit, the interlayer coupling intro¬ 
duces an off-diagonal term to mix the two layers and 
leads to splitting of bands that come from each layer. 
Certain complications arise from the difference in crys¬ 


tal orientation of the two monolayer units. The k • p 
hamiltonian for the bilayer unit at T reads: 

% r '“ H (fc) = 92 + a 2 |&|”(<?3 + 94,tfx) + (23) 

The interlayer coupling mixes the states on two layers 
by the g$ term. The coupling also introduces a renormal¬ 
ized constant energy shift, g 2 , for the central position of 
these bands, instead of go in a single layer unit. This shift 
comes from the virtual hopping process to other bands 
which are integrated out. In this case, this is mostly from 
the hybridization between the top valence bands at T and 
number four band below on the other layer with po char¬ 
acter (see Table IV). The energy shift is positive from 
second order perturbation theory and is proportional to 
the square of the overlap matrix element and the inverse 
of the energy difference. The effective masses are renor¬ 
malized as well in the presence of interlayer coupling and 
these effects are incorporated in the values of the param¬ 
eters go and g 4 . 

At the K point, the orbital character and k • p forms 
for these two layers are different and these can be ac¬ 
counted for by an xz mirror operation on the basis and 
the hamiltonian. Without the interlayer coupling, the 
hamiltonian for a spinless 2H bilayer is: 

?4 v2H (fc) = y (1 + a z ) + fia(rk x a x + k y fi z a y ) 

+ a 2 |fc| 2 (/ 2 + f 3 a z ) + / 4 a 2 (fc 2 d- x - k 2 y cr x - 2rk x k v jl z ay) 

(24) 

The form for spin-orbit coupling is: 

a ^k, 2H = / 5 r/t z ( 1 ^ Z )s z + f 6 Tp z ( 1+ 2 az )sz (25) 
and the interlayer coupling term is: 

A n K ’ m = / 7 ( 1 ^)A, (26) 

We note that in spinless form, only the valence bands 
will split while the conduction bands remain doubly de¬ 
generate by symmetry. The numerical values of the co¬ 
efficients for these k • p hamiltonians are given in Table 
VI, including GW results for spinless monolayer MoS 2 . 


VII. APPLICATIONS 

In what concerns applications of the hamiltonians we 
have derived here, optical transitions 13,33 and Berry cur¬ 
vature effects 48 are two important factors that control 
low-energy dynamics of electrons at different valleys. To 
incorporate interband coupling effects, we use the spin¬ 
less eleven-band TBH to evaluate these physical quanti¬ 
ties. The generalization to the spin case is straightfor- 
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ward. 


A. Interband Optical Transition 

Light couples to the Bloch bands through the gauge 
field. A heuristic way to derive the optical transi¬ 
tion matrix is to start from the expansion of the tight- 
binding hamiltonian at an arbitrary fco, H(k + fco) ss 
H(ko) + k ■ dH(k)/dk\k 0 . The gauge field A enters in 
the replacement k — ► k + eA/h with the electron charge 
— e. The interaction term with the gauge field is V(k) = 
A x V x (k) + AyVy(k) where 7A(fc) = ( e/H)dH(k)/dk l 33 
and A is determined by the polarization of light. This 
is the current operator for the hopping terms. For a± 
circularly polarized light, V±{k) = V x {k) ± iP y (k). The 
interbancl optical transition rate, obtained from Fermi’s 
golden rule, between valence and conduction bands is de¬ 
termined by the matrix element of V(k) between these 
two Bloch wavefunctions ifa: 

J(w) ~ w -2 ^2 | < V’cfcl'PWIV’vfc > 1 2 5(e c (k)-e v (k)-huj) 

k,c,v 

(27) 

The crystal momentum k is the same for these two Bloch 
states in the optical limit and the u>~ 2 factor is from the 
conversion between A 2 and the incoming light intensity 

E 2 . 

In Fig. 9(a), we compute I+(ui) numerically for the 
integrated <t + absorption over the entire Brillouin zone 
and compare with the joint density of states (JDOS), 
Efc c v &{ e c{k) — e„(fc) — hui), between conduction and va¬ 
lence bands. Light absorption starts at the band edge 
between the top valence band and lowest conduction 
band. Because of the broken inversion symmetry, the 
two valleys at K± respond differently to circularly polar¬ 
ized light 49,50 . At K + the orbital character of the highest 
valence band and the lowest conduction band is c ?2 and 
do respectively. a± polarized light adds angular momen¬ 
tum 8m z = ±1. Under 77-3 symmetry, only er + can have 
a non-zero matrix element for the optical transition due 
to angular momentum conservation. To demonstrate the 
optical chiral selection rule in our hamiltonian, we com¬ 
pute the circular dichroism r](k) i9 , given by: 


_ m i, (fc)i 2 -i^-(fc)i 2 

V[ ’ \Vl v {k )\ 2 + \Vl v (k )\ 2 


(28) 


This is the fc-resolved quantity to distinguish cr + and cr_ 
light absorption between the highest valence and the low¬ 
est conduction bands with the matrix element given by 
V±{k) =< ^ c k \P±{k)\ipl >• The photon energy is im¬ 
plicitly chosen to match the energy difference between 
the two bands at each k point. In the inset of Fig. 9(a) 
we plot rj(k ): at K-t, it is either +1 or —1 due to the 
optical chiral selection rule. Near the T point, the lowest 



FIG. 9. M 0 S 2 TBH: (a) <t+ absorption and comparison with 
joint density of states (JDOS);the inset shows the fc-resolved 
degree of polarization jj(fe). (b) Berry curvature in units of 

A 2 . 


conduction band is an odd state and the optical transi¬ 
tion from the highest valence band vanishes due to the 
symmetry selection rule. However, we find a strong de¬ 
gree of optical polarization over a large region near K± 
where the selection rule is not obeyed. In reality, the de¬ 
gree of polarization is degraded by inter-valley scattering 
which is beyond the scope of the present discussion. 


B. Berry Curvature 


In the semi-classical equation of motion for a wave 
packet, transport properties are controlled by the band 
dispersion and the anomalous velocity contribution from 
the Berry curvature fl ra (fc ) 18,48,51 defined as: 


^n(k) — iz • (77fc^nfe) ^ (77 k^nk) 

0 /^\2 \ ' ^ Unk I'Px (k) | Vj n i\z X U n 'k\^y {k) > 

2 ^ ( £ n ( fc )- en ,( fc)) 2 

(29) 

The integral of the Berry curvature over the entire 2D 
Brillouin Zone manifold gives the integer Chern num¬ 
ber, also known as the TKNN number 52 , for that band. 

In Quantum Hall systems the sum of these integers for 
the filled bands corresponds to the Hall conductivity. 
Though the integral is zero for a time-reversal invari¬ 
ant system, the Berry curvature in TMDC monolayers 
is non-zero at generic k points due to broken inver¬ 
sion symmetry. Time-reversal symmetry dictates that 
Qn(—k) = —f l n {k). Based on the TBH model we de¬ 
rived, the Berry curvature can be readily evaluated 53 . In 
Fig. 9(b) we plot the Berry curvature for the top valence 
band of M 0 S 2 . The curvatures at the K± valleys are op¬ 
posite to each other which means that the charge carriers 
in the two valleys will drift differently under an applied 
external electric field. 
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VIII. CONCLUSION 


Appendix A: TMDC TB-Model Simplifications and 
Numerical Parameters 


In this work, we derived ab-initio tight-binding harnil- 
tonians for TMDCs based on a Wannier transforma¬ 
tion. The model hamiltonians retain faithfully the over¬ 
lap matrix elements, the orbital information and the ac¬ 
curacy of DFT calculations on which they are based. The 
eleven-band TBH with first-neighbor and partial second- 
neighbor couplings captures the features of the hybridiza¬ 
tion of the atomic p and d orbitals and reflects the A4i, 
A4 2 and 723 symmetries in a monolayer unit. The spin- 
orbit coupling enters as the on-site LS terms related to 
individual atoms. 

For multiple stacked TMDC layers, we determine the 
interlayer couplings between p orbitals on adjacent X lay¬ 
ers. The two-center Slater-Koster approximation works 
well for modeling interlayer coupling, with hopping terms 
Vp Pi(T (r) and V pP)7r (r) expressed through a simple yet 
transferable empirical function of the pair distance. This 
interlayer interaction term is an essential ingredient in 
understanding TMDC heterostructures when there is 
a twist in the relative orientations of adjacent layers 
or translations between monolayer units. It would be 
interesting to generalize this interaction to other two- 
dimensional layered materials such as graphene and hBN 
where the pair orientation dependence is needed, in ad¬ 
dition to the pair distance. 

Another way to utilize the FTBH is to perform a k 
• p expansion for the low-energy bands at T and Kj-. 
The second-order k ■ p hamiltonian respects the crystal 
symmetry and is constructed to investigate the effects of 
spin-orbit coupling and interlayer hopping in the bilayer 
unit. 

As examples of applications of our TBH, we inves¬ 
tigate the optical absorption and the Berry curvature 
of M 0 S 2 . The tight-binding hamiltonians can form the 
basis for further theoretical investigations, many-body 
physics, and simulations for potential applications under 
external electric or magnetic fields in finite-size nano¬ 
structures 54,55 , in either monolayer or heterostructure 
forms. 
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The 86 parameters for all a , t ■ , t \j , t ■ j and tf'j 

can be simplified to 36 independent parameters, which 
include a subset of e*, and all t\ 5 j , by symmetry 
considerations. Under M 2 and 72-3 symmetries, the hop¬ 
pings to the symmetrical positions are not mutually in¬ 
dependent. Here we summarize the simplifications of the 
truncated tight-binding parameters and tabulate their 
numerical values in Table VII and VIII. For the sets 
of indexes: (a = l,/3 = 2), (a = 4,/3 = 5,7 = 3), 
(a = 7, /3 = 8,7 = 6 ), (a = 10, /3 = 11,7 = 9) with 
the first superscript index corresponding to (+) and the 
second to (—), we have the following relations: 


^■OL t/3 

+(2) = Id 1) I 3 +(!) 

u a,a 4 / 0 .“ ' 4^(3./5 

A 2 ) _ 3 di) , Ido 

l f>,P - 4 r a,« + 4 Z 73,/3 

d 2) = d 1) 

7,7 7,7 

(2,3) = ± ^ 3 dl) _ 1 dD 
, '7,/3 - 1 - 2 7>« 

d2,3) _ , VS/dl) ,(1) 3 dl) 
6 a ,/3 ~ ^ 4 


a,(3 


(Al) 


d 2,3) = Idl) , AjL(1) 

“'7,a 2 ’ ,a 2 

while for (a = 1 ,[3 = 2, a' = 4, f3 r = 5 , 7 ' = 3), (a = 
7,/3 = 8 , a! = 10, j3' = 11, 7 ' = 9), we have: 


dR _ 1(5) 3 ( 5 ) 

^ ol' ,a: 

d4) _ 3 (5) 1(5) 

V,/3 “ 4 7,a + 4 V,/3 

d4) _,(4) _ _ A(6) , \/3.(5) 

/3',a ^ L a',ot ' ^ /3',/3 

,(4) _ _\/3 .(5) 

2 7 d 

d 4 ) _ 7(5) 

S' 2’’^ 


R4)_.(5) ,(4) _ — \/3 (5) ,(4) _ -4, 
6 9,6 — 6 10,6 — a 6 11,6> 6 11,6 — o 1 


-1 


(A2) 


°9,6 


^9,6’ UO ,6 


_ (5) 

2 ” 11 , 6 > ‘' 11,6 — 2 t ' 11 ) 6 


Appendix B: TMDC GW Quasi-Particle Calculation 

We performed a GW calculation using the Quantum 
ESPRESSO 56 and BerkeleyGW 57,58 code which corrects 
the band-gap by computing quasiparticle energies to ob¬ 
tain an improved model hamiltonian. The interactions 
between the ionic cores and electrons were modeled with 
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TABLE VII. Tight-binding independent parameters in units 
of eV for M 0 S 2 , MoSe 2 , WS 2 , WSe 2 based on the DFT results. 



MoS 2 

MoSe 2 

ws 2 

WSe 2 

ei = e 2 

1.0688 

0.7819 

1.3754 

1.0349 

£3 

—0.7755 

-0.6567 

-1.1278 

-0.9573 

£4 = £5 

-1.2902 

-1.1726 

-1.5534 

-1.3937 

ee 

-0.1380 

-0.2297 

-0.0393 

-0.1667 

e 7 = e 8 

0.0874 

0.0149 

0.1984 

0.0984 

£9 

-2.8949 

-2.9015 

-3.3706 

-3.3642 

£l 0 = £ll 

-1.9065 

-1.7806 

-2.3461 

-2.1820 

C l,l 

-0.2069 

-0.1460 

- 0.2011 

-0.1395 

°2,2 

0.0323 

0.0177 

0.0263 

0.0129 

~~t m 

b 3,3 

-0.1739 

- 0.2112 

-0.1749 

-0.2171 

~ t ft 

6 4,4 

0.8651 

0.9638 

0.8726 

0.9763 

t {1) 

6 5,5 

-0.1872 

-0.1724 

-0.2187 

-0.1985 

~p) 

z 6,6 

-0.2979 

-0.2636 

-0.3716 

-0.3330 

A 1 ) 

L 7,7 

0.2747 

0.2505 

0.3537 

0.3190 

t {1) 

l 8,8 

-0.5581 

-0.4734 

-0.6892 

-0.5837 

L 9,9 

-0.1916 

-0.2166 

- 0.2112 

-0.2399 

‘' 10,10 

0.9122 

0.9911 

0.9673 

1.0470 

t {1) 

‘' 11,11 

0.0059 

-0.0036 

0.0143 

0.0029 

z 3,5 

-0.0679 

-0.0735 

-0.0818 

-0.0912 

z 6,8 

0.4096 

0.3520 

0.4896 

0.4233 

t (1) 

‘'9,11 

0.0075 

0.0047 

-0.0315 

-0.0377 

l 1,2 

-0.2562 

-0.1912 

-0.3106 

-0.2321 

A 1 ) 

l 3,4 

-0.0995 

-0.0755 

-0.1105 

-0.0797 

A 1 ) 

l 4,5 

-0.0705 

-0.0680 

-0.0989 

-0.0920 

t 6,7 

-0.1145 

-0.0960 

-0.1467 

-0.1250 

t ^ 

Z 7,8 

-0.2487 

- 0.2012 

-0.3030 

-0.2456 

z 9,10 

0.1063 

0.1216 

0.1645 

0.1857 

t {1) 

‘' 10,11 

-0.0385 

-0.0394 

-0.1018 

-0.1027 

AV 

l 4,l 

-0.7883 

-0.6946 

-0.8855 

-0.7744 

^(5) 

l 3,2 

-1.3790 

-1.3258 

-1.4376 

-1.4014 

AV 

l 5,2 

2.1584 

1.9415 

2.3121 

2.0858 

AS) 

l 9,6 

-0.8836 

-0.7720 

-1.0130 

-0.8998 

AS) 

6 11,6 

-0.9402 

-0.8738 

-0.9878 

-0.9044 

Ho ,7 

1.4114 

1.2677 

1.5629 

1.4030 

A5J 

z 9,8 

-0.9535 

-0.8578 

-0.9491 

-0.8548 

AS) 

6 11,8 

0.6517 

0.5545 

0.6718 

0.5711 

,( 6 ) 

^9,6 

-0.0686 

-0.0691 

-0.0659 

-0.0676 

AS) 

‘'ll ,6 

-0.1498 

-0.1553 

-0.1533 

-0.1608 

- 

•—' 05 

-0.2205 

-0.2227 

-0.2618 

-0.2618 

AS) 

‘'ll ,8 

-0.2451 

-0.2154 

-0.2736 

-0.2424 


normconserving pseudopotentials 59 . The exchange cor- 


TABLE VIII. Atomic spin-orbit coupling strength in TBH in 
units of eV/h 2 . 



Mo 

W 

S 

Se 

7 M/X 

A so 

0.0836 

0.2874 

0.0556 

0.2470 



FIG. 10. Rescaled TBH hamiltonian (red lines) and its com¬ 
parison with GW (blue circles) band structure results for a 
M 0 S 2 single layer unit. 

relation energy between the electrons was treated within 
the Generalized Gradient Approximation (GGA) param¬ 
eterized by Perdew, Burke and Ernzerhof (PBE) 31 . The 
semi-core 4 d, 4 p, and 4s states of Mo were taken as va¬ 
lence states for our DFT and GW calculations. All the 
integrations over the BZ were carried out over a 30 x 
30 x 1 uniform mesh of Appoints, and an energy cutoff of 
1900 eV was used to truncate the plane wave basis used in 
representing Kohn-Sham wave functions. The self energy 
(cr) and the dielectric matrix (e^Q^q)) were calculated 
on a uniform mesh of 12 x 12 x 1 fc-points in the BZ 
with a truncated Coulomb interaction. We employed an 
energy cutoff of 272 eV for the dielectric matrix, and 1750 
bands to calculate e^Q^q) and cr. The conduction and 
valence band quasiparticle energies are converged within 
5 meV using the above parameters as discussed by Mal¬ 
one et al. eo . 

The MoS 2 TBH based on GW can be constructed as in 
the DFT case. GW bands have similar shapes as those in 
DFT and a larger band gap between the groups of valence 
and conduction bands. This effect can be incorporated 
mostly by enlarged coupling constants between orbitals. 
Instead of tabulating a new column for GW-TBH pa¬ 
rameters, we introduce atomic on-site energy shifts and 
scaling factors for the M 0 S 2 TBH parameters. Hopping 
terms in TBH are divided into four groups as different 
bonds in Fig. 2. Each type of them are scaled by differ¬ 
ent scaling factors z, tfj^ = zt^J T , tabulated in Table 
IX under the condition to have the right GW band-gap, 
and the right values for the highest valence band energy 
at P. The GW-TBH is compared with the GW calcu- 
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lation in Fig. 10, showing good agreement between the ^ , „ , T 

, TABLE IX. On-site energy shifts m eV and scaling factors for 

tW °' M 0 S 2 GW-TBH based on DFT-TBH. 



Ae M 

Ae x 

1st M-M 

1st X-X 

1st X-M 

2nd X-M 

M 0 S 2 

0.3624 

-0.2512 

1.4209 

1.1738 

1.0773 

1.1871 
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